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Abstract. The Schrodinger equation defines the dynamics of quantum particles which 
has been an area of unabated interest in physics. We demonstrate how simple transfor- 
mations of the Schrodinger equation leads to a coupled linear system, whereby each 
diagonal block is a high frequency Helmholtz problem. Based on this model, we de- 
rive indefinite Helmholtz model problems with strongly varying wavenumbers. We 
employ the iterative approach for their solution. In particular, we develop a precon- 
ditioner that has its spectrum restricted to a quadrant (of the complex plane) thereby 
making it easily invertible by multigrid methods with standard components. This 
multigrid preconditioner is used in conjuction with suitable Krylov-subspace methods 
for solving the indefinite Helmholtz model problems. The aim of this study is to report 
the feasbility of this preconditioner for the model problems. We compare this idea with 
the other prevalent preconditioning ideas, and discuss its merits. Results of numerical 
experiments are presented, which complement the proposed ideas, and show that this 
preconditioner may be used in an automatic setting. 
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grid, Complex-shifted Laplacian (CSL), Complex-scaled Grid (CSG), Quadrant-definite (QD) 



1 Introduction 

Acoustic, electromagnetic or seismic waves can all be modeled by a Helmholtz equation 
with a wave number that has properties specific to the problem area. In some acous- 
tic scattering applications, for example, the wave number is space independent but the 
boundary of the domain can be very complicated depending on the shape of the object. 
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In electromagnetic scattering there are jumps in the material parameters which lead to 
a piecewise constant wavenumber. In a similar way, the wavenumber in seismic waves 
will contain information about the geological layers in the earths crust. Each of these 
problems pose different challenges to the numerical methods. 

In this article, we focus on the iterative solution of the Helmholtz equations with 
a wave number that is specific to models for breakup problems in chemical systems. 
These breakup dynamics are described by a Schrodinger equation that reduces, in the 
energy range of breakup problems, to Helmholtz equation with a wavenumber that is 
continuous in the space variables and can become very large near the boundary of the 
domain. One example is the disintegration into four charged particles of the H2 molecule 
when it is hit with a single photon [1]. 

The prevalent practice for solving this type of problem requires massively parallel 
computers [ ] and they use a significant portion of the resources of large computer fa- 
cilities. The long term aim of this research is to replace this practice by efficient iterative 
methods. 

The Helmholtz equation has often outgoing wave boundary conditions. Fixing ho- 
mogeneous Dirichlet boundary conditions, on the boundaries of the truncated numerical 
domain, leads to artificial reflections in the domain of interest. These reflections are nu- 
merical errors and must be diminished by damping the outgoing waves at the domain 
boundaries. To bring this about in our numerical solution method, we employ exterior 
complex scaled [ ] absorbing boundary layers (henceforth ECS-ABL). There is a long his- 
tory of this type of absorbing boundary condition for chemical reactions [4]. This treat- 
ment is equivalent to the use of perfectly matched layers (PML) [ , ] and leads to a non- 
Hermitian discrete problem [7]. For a review on transparant and absorbing boundary 
conditions for the Schrodinger equaton we refer to [8]. 

For Krylov-subspace methods, the main challenge is to find a good preconditioner. 
Over the years there have been different approaches to preconditioning the indefinite 
Helmholtz equation. One line of research is based on a shifted Laplacian preconditioner 
that started with the work [9,10] ( Bayliss, Goldstein and Turkel). They used the Laplacian 
and the positively shifted Laplacians as preconditioner. 

This was later successfully generalized into a robust method, known as the complex 
shifted laplacian (CSL) preconditioner, by Erlangga, Vuik and Oosterlee using complex val- 
ued shifts in [11, 12]. Introducing a complex shift pushes the spectrum of the Helmholtz 
operator into a region that is favorable for multigrid methods [13-15] to approximately 
invert the preconditioning problem. It is well-known that multigrid efficiency can read- 
ily be exploited only for problems having (positive or negative) definite spectra. In the 
indefinite case [ ], both vital components of multigrid, i.e., smoothing, and coarse grid 
correction suffer severe degradation, and consequently this results in divergence of the 
method [16]. 

An alternative preconditioner that, in addition to shift, also scales the Laplacian was 
derived from frequency shift time integration by Meerbergen and Coyette [ ]. By appro- 
priately choosing the shift and the scale it is possible to restrict the spectrum of precondi- 



3 



tioning matrix into one quadrant of the complex plane. We call this type of preconditioner 
a quadrant definite (QD) preconditioner. 

In [7], we proposed the complex-scaled grid (CSG) preconditioner, and demonstrated 
its utility in connection with indefinite Helmholtz problems constructed with ECS-ABL. 
Both the CSL and the CSG preconditioners have similar performance and are based on 
similar ideas. The CSL translates the spectrum, while the CSG rotates it, thereby placing 
it in a region which is multigrid favorable. Both of these preconditioners depend on the 
translation magnitude or the rotation angle which has to be tuned for specific problems. 

This paper studies a preconditioner based on a scaled translation of the spectrum 
that restricts it to one quadrant of the complex plane. We evaluate it on a set of model 
problems representative for breakup problems that are derived in the paper. While its 
efficiency is found to be between that of the Laplacian preconditioner and the CSL/ CSG 
preconditioners, the main merit is its ease of invertibility by multigrid methods that use 
well-known standard components. This is a clear advantage of using the QD method, as 
for the CSL/ CSG preconditioners, multigrid has to be tuned for different wavenumbers. 
Moreover, a shift for the CSL preconditioner (or equivalently a rotation angle for the CSG 
preconditioner) is apparantly only available through a hit-and-trial rule. In comparison, 
the QD preconditioner may be used in an automatic setting. 

Both the discretisation and the absorbing boundary conditions used in this paper are 
of low order of accuracy. Both can be replaced by higher order methods, however, the 
focus of the paper is on the working of the iterative methods and this can be studied with 
the low order methods since the higher order discretisation and boundary conditions 
have similar spectral properties. 

In Section 2.1, we give the transformation of the Schrodinger equation to a coupled 
Helmholtz problem, and derive the model problems for this study. The details of ECS- 
ABL and the discretization are given in Section 3. Also reviewed here, are the spectral 
properties of the discrete operator. Next, in Section 4, we describe the QD preconditioner 
in detail, and give the multigrid algorithm which we use for approximate inversion of the 
preconditioners. This is followed by numerical experiments, which are given in Section 
5. Some conclusions mark the end of the paper in Section 6. 

2 From the Schrodinger equation to a coupled Helmholtz prob- 
lem 

2.1 The model problem 

In this section we derive the model Helmholtz problem that we use in this paper to bench- 
mark iterative solvers. The model problem is 

' (-A hih -k 2 (x,y))u{x,y)=f(x,y) on [0,fl] 2 CR 2 , 
< u(x,Q)=0 Vxg[0,a] and «(0,y)=0 Vye[0,a], (2.1) 
ABC on u(x,a) Vxe[0,a] and u{a,y) Vye[0,a], 
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where ABC denotes outgoing wave boundary conditions, see Section 3, and 

A W2 =d^+d yy -2 -J (2.2) 

denotes the radial part of the Laplacian in spherical coordinates with l\, I2E N. The 
wavenumber k 2 (x,y) =2m(E — V(x,y)) depends on a potential V(x,y) that varies contin- 
uously in x and y in the domain [0,a] 2 , E > is the energy and m > is the mass of the 
system. The right hand side f(x,y) is assumed to be zero outside [0,b] 2 with b < a so that 
the Helmholtz problem becomes a homogeneous problem in a strip near the boundaries 
with the ABC. 

2.2 The Schrodinger equation 

To derive this model we start from the driven Schrodinger equation 

CH-E)xp(r 1 ,r 2 )= ( p(r 1 ,r2), (2.3) 
with r\, t2 G K 3 and where H denotes the Hamiltonian and is given by 

n = -^ n ^ 1 -^r 2 + Vi{\ri\) + V 1 (\r 2 \) + V ll {r l ,r 2 ) (2.4) 

with V\ and Vj local potentials that only depend on magnitude of r ; . The potential V\2 
depends, usually, on the relative distance between r\ and r%. The mass m > scales the 
Laplacians. The right hand side of (2.3), $(11,12), is assumed zero if |ri | > b or \i2 \ > b, and 
can model an incoming electron that impacts in the system [ ], or alternately, represents 
the dipole operator working on a groundstate if the model is used to compute photo- 
ionization [1]. 

For these breakup problems the solution $(v\,V2) is an outgoing wave in any direction 
similar to the Sommerfeld radiation condition. This leads to a six dimensional problem 
on an unbounded domain. The problem can also be interpreted as a 6D Helmholtz prob- 
lem 

(-^ev-k 1 (r l ,r l )) ip(t lt t 2 ) = f(ri,r 2 ), (2.5) 

where \(r\,ri) = y / 2m(E — V) with V denotes the sum of all potentials. This becomes 
a Helmholtz problem with a constant wave number, k = \/2mE, in the regions of space 
where the potentials go to zero. This 6D problem is hard to solve with the current gener- 
ation of numerical methods. 

2.3 Expansion of the solution in partial waves. 

In this section we discuss the reduction of the 6D problem to a coupled set of 2D prob- 
lems. At large distances the solution behaves as a spherical wave emerging from the cen- 
ter of mass of the system. It is therefore common practice [19,20] to rewrite the equation 
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(2.3) in spherical coordinates. The Laplacian operator then splits into a radial operator 
and the angular operator differential [ ]. The coordinates are written as r\ = (jOi,Qi) and 
?2 = {pi/^i)/ where Q denotes (0,cp). The solution is then written as a series 

00/^00/2 

lp{n>n)=Yj E E E J / ? /im 1 ,Z 2 m 2 (|0l/|02)Y/ imi (Oi)y Zimi (O2), (2.6) 
li=0 mi=-li l 2 =Qm 2 ——l 2 

where Y/ m (Q) are the spherical harmonics, the eigenfunctions of the angular differential 
operator of the Laplacian in spherical coordinates [ ]. In physics this decomposition 
is referred to as the partial wave expansion and the functions ipi 1 mi,i 2 m 2 {pi, pi) are called 
partial waves. 

When this proposal, (2.6), is substituted in (2.3) we find an equation for xp^ mi ,l 2 m 2 (i°i/i°2) 
for all l\ > 0,/ 2 > 0, \mi \ < l\, \m 2 \ < Z 2 that is coupled to all other partial waves. 

1 . 



CO l\ 



V 2 (2.7) 

E E E E V hm 1 hm 2 ^[m[X 2 m' 2 {pl'P2)^l' 1 m' v V 2 m' 2 {Pl'Pl) = <Phmi,km 2 t 



l' 1 =0m 1 =-l 1 l 2 =0m' 2 =-l' 2 

where the coupling potentials are calculated as 

VimnhmrW^ipvpi) =J d0.idC\ 2 Yf im (Qi)Y^ m2 (Q 2 ) 



x[Vi(|r 1 |) + y 2 (|i- 2 |) + y 12 (r 1 ,r 2 )]Y /;m; (Q 1 )Y /2m2 (n 2 ) 



(2.8) 



and (pi imii i im2 is partial wave of the right hand side. 

When the potentials V\, V 2 and V\2 are spherically symmetric the system decouples. 
When it is cylindrically symmetric the different m\ and m 2 are decoupled. But in general 
the system is fully coupled. Furthermore, it is common practice to truncate the infinite 
series in / at a finite l max so that it becomes a finite system of coupled partial differential 
equations. 

The boundary conditions for Equation (2.3) translate in spherical coordinates into ho- 
mogeneous Dirichlet xp(pi,0) =0 for all pi and i/?(0,|0 2 ) =0 for all p 2 . This is typical for ra- 
dial problem since rhOi =0 and p 2 =0 is now the origin of the coordinate system [ ]. The 
outgoing boundary conditions translate then into outgoing boundary conditions pi — > 00 
or p2 — > 00. We will elaborate on this topic in Section 3. 

The partial wave expansion can also be written down for a single particle Hamilto- 
nian. It then involves an expansion over a single angular function Y/ m and subsequently 
leads to a coupled system of ordinary differential equations. On the other hand, the 
Hamiltonians currently studied in the physics and chemistry communities involve three 
or more particles. For three particles the driven Schrodinger equation is a 9-dimensional 
equation that, after expansion in partial waves, becomes a set of coupled 3D PDEs. 
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2.4 Blocked structure and Iterative methods 

The system (2.7) has a very particular structure. Since the differential operators are block 
diagonal in the spherical expansion, they only appear on the diagonal blocks of the equa- 
tion. The blocks are only coupled by the potentials defined in equation (2.8). The Hamil- 
tonian W can be written in blocked matrix notation as 

vi imi i 2 m2;hmilim-i ~ 2m hth~^~ *hm\h.mi}kmihm2~ ^ ••• , (2.9) 

where the A; / 2 are the radial differential operators defined in (2.2). This can be written 
as a coupled Helmholtz operator 

- k hm 1 l 2 m 2 ;hm 1 l 1 m 1 (P^P2) ~ A h,h -\m x \im^m x hmt (Pl'Pz) •■• I (2.10) 



After discretization of the differential operators on a grid discussed in detail in Section 
3, we arrive at a system of linear equations, Ax = b. The matrix A will have the same 
blocked structure as the coupled system of partial differential equations above and we 
can write: 



(A n 


A U A13 


...\ 


fxi\ 




An 


A22 






b 2 


A 2 i 


A33 




*3 


= h 


V- 




■••/ 


{'■/ 


v) 



(2.11) 



where the discretized differential operators will only appear in the diagonal blocks An. 
Since the differential operators will lead to the largest eigenvalues, the condition number 
of the full matrix A will also be determined predominantly by the diagonal blocks An. 
After discretization of p\ on n grid points and P2 on n grid points, a single block is a 
sparse matrix of size n 2 xn 2 . 

The solution method for solving this type of breakup problems as employed in [1] 
is iterative. This method was developed in [22] and exploits the particular block struc- 
ture of A. A block diagonal preconditioning matrix M is constructed that contains only 
the diagonal blocks A,-,-. Since the largest eigenvalues and eigenvectors of M and A are 
very similar, M _1 A has a smaller condition number, and therefore, M _1 proves to be 
a good preconditioner for any suitable Krylov-subspace method. However, note here 
in particular, that the strategy in used in [1] is to exactly invert the blocks (each of size 
n 2 x n 2 ) within the preconditioning step. Inasmuch as each diagonal block represents a 
two-dimensional system, the diagonal block matrices can be inverted possibly on a single 
processor. The coupled system, however, requires the inversion of many diagonal blocks 
and requires a cluster. 
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However, the problems currently under investigation in the physics and chemistry 
communities such as the impact-ionization problems or problems where electronic mo- 
tion and nuclear motion are combined described in section 2.5.3, each diagonal block 
constitutes a three dimensional problem and renders itself too unwieldy for exact inver- 
sion. 

We therefore study in this paper the multigrid-preconditioned iterative solution of 
the diagonal block only, and not the entire problem as a whole. The diagonal blocks (2.9) 
corresponds closely to the model problem (2.1). It is important to understand that the 
complete process now involves two independent iterative schemes, the outer scheme for 
approximately inverting the entire system via a preconditioned Krylov process, and the 
inner iterative scheme which uses multigrid preconditioning for approximately inverting 
the diagonal blocks within the outer preconditioner. The latter alone forms the subject 
matter of this paper. We have chosen to restrict the dimensions to two for this study. 



2.5 Examples 

To illustrate the significance of the coupled system of partial differential equations we 
give a few example physical systems that are currently studied with the approach. We 
cite the relevant papers. 

2.5.1 The dynamics of two electrons in a Helium atom 

The Helium atom is a quantum system that has two electrons with a negative charge and 
one nucleus which has a positive charge of unit two. Since the nucleus is much heavier 
than each electron, the position of the nucleus is taken as the center of the coordinate 
system. In this coordinate system the first electron is at r\ and the second electron at r 2 
The potentials in equation (2.3) are then 

^lOi) = -A, V 2 (r 2 ) = -^, V 11 {r lr r 1 ) = v ^—^ (2.12) 
\n\ \r 2 \ \n-r 2 \ 

To arrive at the potentials in the coupled problem (2.8) the multipole expansion 

1 -£4i^(c°s(M) (2-13) 



\n-r 2 \ yp 



> 



is used to expand V\ 2 . Where p< and p> denote, respectively, the smallest and largest 
of p\ and p 2 . The angle 6\ 2 is between the vectors r\ and r 2 . Since V\ and V 2 are central 
potentials they will appear as —2/p\ and —2/p 2 on the diagonal blocks of (2.8) when 
l\ = l[ and h = l[. The multipole expansion (2.13), however, will lead to potentials that 
couple the blocks with different / values in equations (2.7). Since the problem is sym- 
metric around the z axis, different m blocks are decoupled. Recent processes in Helium 
studied with this approach are one and two-photon double ionization [23,24]. 
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2.5.2 The dynamics of two electrons in the Hydrogen molecule 

The Hydrogen molecule consist of two negatively charged electrons and two protons 
with a positive charge. The two protons are much heavier than the electrons. After 
the Born-Oppenheimer approximation the two protons can be considered fixed in space. 
The dynamics of the two electrons are governed by equation (2.3), where the potential is 
given by the static field of the charged protons. If we take a coordinate system around 
the center of mass of the protons and R is the vector connecting the two protons and r\ 
and r 2 the coordinates of the electrons, the potentials in (2.3) are 

yi(ri) = -^W-R^727' 72( * )= -^^-mW (2 - 14) 

V 12 (r lf r 2 )= . , (2.15) 
\ri-r 2 \ 

The first is the attraction of the first electron to the two protons. The second is the same at- 
traction but for the second electron. The third potential is the electron-electron repulsion 
because both have a negative charge. To derive the potentials in the coupled basis we 
use, again, the multipole expansion (2.13) for each of the potentials. Now all potentials 
couple the blocks. Again, an example of a process studied in this approach is one-photon 
double ionization [1]. 

2.5.3 Impact ionization of Helium or the Hydrogen molecule 

When an additional electron with sufficient energy collides and breaks up the Helium 
atom or Hydrogen molecule that have already two electrons, we are tracking three par- 
ticles. We then have a 9D problem. If we denote the coordinate of the third, impacting, 
electron as r$ we end up with Helmholtz operator 

- A ri - A r2 - A r3 ~k 2 (ri,r 2 ,r 3 ) . (2.16) 
After the partial wave expansion we end up, again, with a coupled problem 

i a 2 i a 2 i a 2 , h(h+i) , h(k+i) h(h+i) 



2m dpi 2m dpi 2m dpi 2p\ 2p 2 2p\ 

+ V(p 1 ,p 2 ,p 3 )-E ip(pi,p 2 ,p 3 )+ E Q l/2w; ^(pi,p 2 ^3)^ / /^(pi^2^3)=0 

(2.17) 

where V(pi,p 2/ p 3 ) is a diagonal block potential while C{pi,p 2 ,pj,) couples the blocks. Be- 
cause of the scale of these problem there are currently no converged results for this prob- 
lem. A similar high dimensional problem can be formulated for Hydrogen molecule 
when no Born-Oppenheimer approximation is applied. Then the motion of the electrons, 
r\ and r 2 , is coupled to the motion of the protons R. Solving these problems is a great 
interest in the scientific community. 
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3 Discretization 

3.1 Absorbing boundary conditions 

In order to solve equations, such as the Helmholtz equation, defined on an unbounded 
domain Qo C R rf numerically, an assumption is made on the asymptotic behavior of 
the solution. The truncated computational domain is a bounded subset Q C Qo of the 
original one, with artificially introduced boundary conditions that imply the postulated 
asymptotic behavior. A commonly used example are the first order Sommerfeld radia- 
tion boundary conditions applied to the homogeneous Helmholtz problem, jj| +iku = 0, 
where n is the outward normal on the boundary 9Q. An exponential decay of the solution 
depending on the constant wave number k is assumed towards the boundary. 

In more complicated Helmholtz models such as the one derived from the Schrodinger 
equation (2.1) more robust techniques are preferable. In the perfectly matched layer (PML) 
approach [5] a small boundary layer r C R is added beyond any point of domain trun- 
cation. On this finite layer the continuous model is adapted to capture the asymptotic 
behavior, with trivial boundary conditions at the end of the layers dT. This idea is equiv- 
alent to a complex coordinate stretching [6,25,26] in the boundary layers, where the orig- 
inal equation is used in the new coordinate system T z C C d with homogeneous Dirichlet 
boundary conditions at the end 3r z , also known as exterior complex scaling (ECS) [3,27]. 
In general we can define an analytic continuation on the layers by 

, J x, xGQ; 

zw ~\ x+if{x), xer, 

with / S C 2 increasing (e.g. linear, quadratic, . . . ) and lim f(x)=0. We denote the image 

of the layer T z = z(T) and call it the complex contour. This boundary layer method does 
not need an explicit input of the wave number and it can easily be tuned in numerical 
experiments. Because of the straightforward mathematical meaning the ECS method is 
interesting in numerical analysis. 

3.2 Finite difference 

ECS boundary conditions and their application in chemical reactions have been used in 
finite difference, B-spline and spectral element discretization [28]. Finite Elments meth- 
ods are hardly used for this type of problems because the computational domain is often 
a square or a rectangular strip. In this article we use finite differences since this low order 
discretisation can already help us to understand the convergence of the iterative method. 
We define a one-dimensional uniform grid (z/)o</<n on the real interval [0,1] with Zq = 
and z n = 1 and mesh width h = 1/n S R. Starting in 1, we apply linear ECS, so the ab- 
sorbing layer is a line connecting 1 and K Z GC henceforth denoted by [1,R Z ]- A second 
uniform grid (zj) n <j< n+m discretized this complex contour, with z n+m = R z and complex 
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mesh width h 7 = (R z — l)/m. The union of these two grids is the ECS grid 

{zj) <j< n+m on [0,1}U[1,R Z ] (3.1) 
in the entire ECS domain. We will denote the fraction 7 = h 7 /h = (R z — 1 ) / (R — 1 ) . 



0.2 



D) 0.1 

to 

I 





-0.1 



+ z 

+ z n+m 
n+m-1 



+ + + + + + + + + + + + + + + + + + + + + + + + Z 

1 2 n-Tn 



x=0 



x=1 



x=R 



real 



Figure 1: Discretized ECS domain Zj. The ECS domain is discretized with complex mesh widths on the complex 
contour. 



A thorough numerical analysis of the negative Laplace operator L= — A discretized on 
this ECS domain yields some important insights for the use of ECS on more general oper- 
ators. To approximate the second derivative we employ the following standard formula 
for un-equal mesh sizes, and non-uniform grids: 



d 2 u 



(*;)' 



1 



1 



1 



for non-uniform grids in grid point where and hj are the left and right mesh widths 
respectively and may belong either to the h category or to the h 7 category. The formula 
can be easily derived as an exercise in Taylor expansions and it reduces to regular second 
order central differences when h;-\ = h;, i.e., in the interior real region (0,1), and in the 
interior of the complex contour (1,R Z ) because the scaling function / is taken to be linear. 
The only exception is the point z n where at most we lose an order of accuracy, however 
with ample discretization steps, the overall accuracy is anticipated to match up to second 
order. We will denote the resulting discretization matrix L/,. 



3.3 Spectral properties 

The hardest model problem, from an iterative point of view, is the problem with l\ = 
and h = since the problem is then at its most indefinite state. For larger l\ and h the 
problem becomes more definite. Therefore we focus on the remainder of the paper on the 
problem with l\ = and h = 0. 

The spectrum of the discretization matrix L/, determines the convergence behavior of 
iterative methods such as Krylov subspace methods and multigrid schemes for solving 
any system L h Uh = \. The spectrum cr(L h ) is drastically different from the spectrum 
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c(L) of the continuous operator, on the undiscretized ECS grid [0,1]U[1,R Z ]. Indeed, 



2 



c(L) = < I ^- ) \j G Nq f , is an infinite set of points on the complex line pe t2e ", with pGR + 



and Oct the complex angle of the complex boundary R z . The shape of the spectrum cr(L^) 
of the discretization matrix is less obvious as follows from the next lemma, that is proved 
in [7]. 

h 

Lemma 3.1. Consider the ECS grid (3.1) and the discretization matrix L;,. Define 7 = 
Then the eigenvalues of L/, are the solutions of 

F / A x_ ten(2np(A)) cos(p(A)) 
W ~tan(2 mt? (A))^cos(^(A)) ' 

with p(A) = 2arccos(l- f h 2 ), q(A) = iarccos(l- § j 2 h 2 ). 

For the Laplace problem the ECS discretized spectrum has the typical Y-shape of a 
pitchfork. There is a clear complex branch associated to eigenvectors located on the com- 
plex contour, along the complex line [0,4/ 'h 2 ], and a branch closer to the real line [0,4/ 'h 2 ] 
that corresponds to eigenvectors located on the real domain. The smallest eigenvalues, 
in the small tail of the pitchfork, belong to the smoothest eigenvectors spread over the 
entire ECS domain. They lie close to the smallest eigenvalue of the continuous ECS oper- 
ator L (Figure 2), that is along the complex line [0,4/h 2 ], where h a =R z / (n+m + 1) is the 
complex mesh width, belonging to a straight complex grid connecting and R z . 



4 The QD preconditioner and multigrid 

4.1 The preconditioner 

We use a preconditioning operator that has a spectrum bounded by a single quadrant 
such that it can be approximately inverted with standard multigrid components, which 
are clearly unstable for indefinite problems. 

In this article we compare the use of a preconditioner Z which is a scaled and shifted 
version of the original Helmholtz operator Z= —Ai lr i 2 — k 2 (x,y) defined in equation (2.1). 
We propose to use 

Z = 5 2 Z + (l-i5k) (4.1) 

where 5 G IR is chosen such that Z is definite. This preconditioner is very similar to the 
one proposed in [17]. Suppose Ao is that eigenvalue of the original operator Z which 
has the smallest real part. We can then choose 5 such that — S \Aq\ +1 >0. For a Helmholtz 
problem with a constant wave number k and l\ =0 and I2 =0 this would mean that 5<l/k. 

The eigenvalues of the preconditioned operator Z~ X Z lie inside a circle of radius 
1/ (2S 2 ) \ l — (kS) \ centered around i (| — i^)- We can readily see this with the follow- 
ing arguments. The preconditioner Z = S 2 Z + (1 — iSk) has the same eigenvectors as the 
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Figure 2: The eigenvalues of the ECS Laplacian discretization matrix (•) lie along a pitchfork with a Y 
shape, a result in Lemma 3.1. A part of the eigenvalues lie close to the eigenvalues the same Laplace problem 
restricted to the interior real domain (<). In a similar way part of the eigenvalues lie close to the eigenvalues 
when the Laplacian is restricted to the complex contour (>). The inset shows the area around the origin where 
the smallest eigenvalues are approximated by the smallest eigenvalues of the Laplace problem defined on the 
complex line [0,R 2 ] (A), (color online) 



Z. The eigenvalues of the preconditioned system Z X Z are therefore given by 



where A is an eigenvalue of Z. We assume that the eigenvalues of Z are located in the 
lower half of the complex plane, cr{Z) C C . Then a{Z~ l Z) is inside the circle that is 
the image of the real axis of the transform (4.2). This circle goes through 0, l/i5 2 6lR 
and —if (kS 3 ) G zTR, so the center c is the crossing point of the lines 9ft (z) = 1/ (2<S 2 ) and 



And so the radius is r = \c\ = 1 / (25 1 ) \l-if (kS) \ . 
4.2 Multigrid 

Heuristically we note that multigrid (for convergence and efficiency), has more stringent 
requirements on the condition number of the spectrum when it crosses into different 
quadrants of the complex plane, than when it does not. The preconditioner in Section 4.1 
has the property that its spectrum is restricted to the fourth quadrant. It can therefore be 



A 



(4.2) 



S 2 \+(l-iSk)' 



S(z) = -l/(2^ 3 ), 
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Figure 3: An F|(l,l) cycle, obtained by Algorithm 1, with 7^ = 1 and 7 C =4. o stands for smoothing, \ stands 
for restriction to a lower level, / stands for prolongation to a higher level, • stands for exact solution. 



very efficiently inverted by multigrid using the standard components, which include co- 
Red Black Jacobi, with a> = 1.05, Full Weighting averaging for restriction, Bilinear interpo- 
lation for prolongation, and rediscretization on the coarse grids; in a simple V(l,l) cycle 
set-up. For experimental purposes we compare the performance of this quadrant-definite 
(QD) preconditioner with the CSG and the CSL preconditioners. In [7], we showed that 
the CSL and the CSG preconditioner can be inverted efficiently using multigrid based 
on matrix components, such as ILU-smoother and the Galerkin coarse grid operator in a 
V(0,1) cycle set-up. This study is more focused on using matrix-free components, which 
leaves only the o;-Jacobi smoother, and the discretization coarse grid operator for the CSL 
and the CSG preconditioners. Moreover, this multigrid has to be employed in an F^(l,l) 
cycle set-up. 



Algorithm 1 is a slight variation of the standard multigrid algorithm in [ ], and 
yields various cycle types depending on the values of and 7 C in a unified manner. 
In this algorithm, / indicates the current level, C the coarsest level, and A, the discrete op- 
erator (i.e., CSL/ CSG/ QD precond. operators) at various levels, bi is the right hand side, 
uf is the starting guess, and V\ and V2 are the number of pre and post smoothing sweeps. 
jf and 7 C are the cycle indices used at the fine and the coarse levels, respectively. E.g., 
calling this method with 7^ = 1 and 7 C = 1 gives the standard V cycle, 7/ = 2 and 7 C = 2 
gives the standard W cycle, while 7/ = 1 and j c = 2 yields the standard F cycle, jf = 1, 
and 7 C = n renders an F cycle with n — 1 recursions on the coarse levels. We found that 
in the context of inverting CSL and CSG based Helmholtz preconditioners, F cycles with 
2 and 3 recursions on the coarse levels were particularly beneficial over the standard F 
cycle. F cycle with 3 recursions is abbreviated as F^ and is shown in Figure 3. 



14 



Algorithm 1 Multigrid pseudocode 



u1 l+1 =MG(l,Ai,buuf,v 1 ,v 2 ,C^ f , lc ). 

(0) Initialization 

- if / = C, = exact (A;,fy); Bail out ; endif 

- Build the coarse-grid operator A/ + i, and the restriction and prolongation 
operators. 

(1) Pre-smoothing 

- Compute wj" +1/3 by applying Vi(>0) smoothing steps to u™: 
u™ +1/3 = smooth 17 ! (Ai,b h u?) . 

(2) Coarse grid correction 

- Compute the residual r" 7 = b/-A/u™ +1/3 . 

- Restrict the residual ff +x = l\ +l ff . 

- Compute the approximate error 

ef +1 from the defect equation. A/ + i e^ 1 =Fj^ 1 

by the following recursion 



if l = C, ef +1 = exact {A 1+1 ,rf +1 ); endif 
If 1<C, approximate ef +1 recursively: 

e™' 1 =0- 

do /' = 1 to 7 C 
if i==l, 

g^ 1 = MG(/ + l / A /+1/ rf; i/ ^ 1/ i/ 1/ t/ 2/ C /7// 7 (; ) 
else 

g^- 1 = MG(/ + l / A ;+1/ r^ 1/ ^ 1/ t/ 1/ t/ 2/ C / 7 C/ 7 / ) 
endif 
continue i 
endif 



- Interpolate the correction ^? = A+\ ^T+i 



Compute the corrected 



approximation on Q, uf +2/3 = uf +1/3 

(3) Post-smoothing 

- Compute by applying V2(>0) smoothing steps to u" t+2/3 : 



uj 1 =smoo 



, th V2( M «+V2 /Al/fel) 
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Preconditioner 


Multigrid 


mg eye. 
per prec. 


Bi-CGSTAB 
iter , cputime 


eye, smooth., co 


mg-conv. , # cycles 


CSL 

06^2 = (-1,-0.3) 


V\(l,l), a;-Jacobi 
a; = 0.8 


0.43 , 17 
4.21s 


1 


60 , 2m lis 


CSG 

e a = T/i4 


¥((1,1), a;-Jacobi 
o; = 0.8 


0.39 , 15 
3.18s 


1 


62 , 2m 2s 


QD 

Re(A ) = -2.6xl0 4 


V(l,l), co-RB Jacobi 
a; = 1.0 


0.09 , 6 
1.2s 


1 


170 , 5m 39s 



Table 1: Multigrid performance and comparison of the three precondtioners for MP1 with fc=160, 256 cells in 
the interior region, and 64 cells in ECS-ABL on all four sides of the domain. ECS angle used is t/6. 



5 Numerical Results 



In this section we conduct numerical experiments on the Helmholtz model problems 
given by the generic prototype Zu(x,y) =f(x,y). Again we focus on problems with l\ =0 
and I2 = because these are the hardest problems. The operator Z is defined in each of 
the following case as follows: 

Z = -A-k 2 MP1 (5.1) 

Z = -A-v(^ + ^)-k 2 MP2 (5.2) 

Z = -A-----k 2 MP3 (5.3) 

x y 

For MP1, f(x,y) is the Dirac delta function that stays zero throughout the domain 
save one point in the middle where it assumes the value 1. For MP2 and MP3, f(x,y) = 
e x2+ y 2 . We use ECS-ABL on all four sides with MP1, and on the north and the east sides 
only with MP2 and MP3. The model problems are solved iteratively with Bi-CGSTAB 
preconditioned with multigrid approximated inverses of the following operators: 

M CS L = -A + {Pi + iP2)k 2 , (5.4) 

M CSG = Z, on the grid rotated by angle 6 a in the complex plane, see [ ] (5.5) 

1 

J\Aqd = (1 — 7— — tt — rrZ Ao is the eigenvalue of Z with the smallest real part. (5.6) 

|Re(AoJ| 

Mqsl is the Complex-shifted Laplacian as appears in [12]. A small complex shift is added 
to the Laplacian operator. This imparts a rectangular translational effect on the operator 
spectrum. Mqsg is the original operator discretized on the so-called Complex-scaled Grid 
and appears in detail in [ ]. The basic mesh size has been multiplied with e l *. This 
imparts a rotation to the operator spectrum about the origin by an angle equal to 6 a . 
Mqsg is as efficient as Mqsl i n general, and slightly better for the current problems. 
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Preconditioner 


Multigrid 


mg eye. 
per prec. 


Bi-CGSTAB 
iter , cputime 


eye, smooth., to 


mg-conv. , # cycles 


CSL 

(01,02 = (-1,-0.4) 


Pf(l,l), <x>-Jacobi 
a; = 0.5 


0.53 , 22 
13.4s 


1 


137 , 7m 34s 


CSG 

Qa = n ln 


F^(l,l), a;-Jacobi 
a; = 0.5 


0.53 , 22 
14.4s 


1 


143 , 7m 36s 


QD 

Re(A ) = -16.88 


¥(1,1), cv-KB Jacobi 
a; = 1.0 


0.15,8 
5.2s 


1 


357 , 19m 40s 



Table 2: Multigrid performance and comparison of the three precondtioners for MP2 with A = 7, fc = 4. The 
domain is a square of 50 units. 512 cells are used in the interior region, and 128 cells in ECS-ABL on the east 
and the north side of the domain. ECS angle used is 




(a) Solution computed in Exp. 1 (b) Solution computed in Exp. 2 




50 10 60 70 ' 

60 80 

(c) Solution computed in Exp. 3 (d) Solution computed in Exp. 4 

Figure 4: All these four solutions were computed for each of the 4 numerical experiments listed in the tables. 
The first solution as spherical waves ensuing out from the domain center. The other three solutions show 
evanescent waves, also known as single ionization, near the west and the south boundaries of the domain. At 
these edges the spatially dependent wavenumber grows exponentially in the model problems, (color online) 
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Preconditioner 


Multigrid 


mg eye. 
per prec. 


Bi-CGSTAB 
iter , cputime 


eye, smooth., to 


mg-conv. , # cycles 


CSL 

(01,02 = (-1,-0.6) 


Pf(l,l), o;-Jacobi 
a; = 0.8 


0.32 , 13 
6.45s 


1 


60 , 3m 9s 


CSG 

a = */l3 


Fj(l,l), a;-Jacobi 
a; = 0.8 


0.32 , 13 
6.3s 


1 


61 , 3m 10s 


QD 

Re(A ) = -4.19 


V(l,l), to-RB Jacobi 
a; = 1.05 


0.17,8 
1.2s 


1 


164 , 9m 


Table 3: Multigrid performance and comparison of the three preconditioners for MP3 with k = 2. The domain 
is a square of 50 units. 512 cells are used in the interior region, and 128 cells in ECS-ABL on the east and the 
north sides of the domain. ECS angle used is n ld. 


Preconditioner 


Multigrid 


mg eye. 
per prec. 


Bi-CGSTAB 
iter , cputime 


eye, smooth., to 


mg-conv. , # cycles 


CSL 

{h>fr = (-1,-0.6) 


a;-Jacobi 
a; = 0.8 


0.32 , 13 
15.8s 


1 


210 , 18m 20s 


CSG 

9a = n : /l3 


F^(l,l), a;-Jacobi 
a; = 0.8 


0.31 , 12 
14.6s 


1 


160 , 14m 14s 


Re(Ao) = -16.18 


V(l,l), a;-RB Jacobi 
a; = 1.05 


0.13,7 
9.4s 


1 


545 , 46m 40s 



Table 4: Multigrid performance and comparison of the three precondtioners for MP3 with fc = 4, The domain is 
a square of 75 units. 768 cells are used in the interior region, and 128 cells in ECS-ABL on the north and the 
east sides of the domain. ECS angle used is 



Numerical experiment results are reported for multigrid invertibility of the precondi- 
tioner and the observed efficiency of preconditioned Bi-CGSTAB. Multigrid invertibility 
is reported as the average multigrid convergence factor (mg-conv.) and the total num- 
ber of cycles that the algorithm required to converge for the preconditioner taken as a 
stand-alone problem, mg-conv. is actually the geometric mean of the observed resid- 
ual decay rates during multigrid cycles, computed over the last 5 cycles. The CPU-time 
is also reported. Bi-CGSTAB efficiency takes into account the number of iterations of 
the algorithm for convergence. Note that each Bi-CGSTAB iteration has two embedded 
multigrid cycles for preconditioning, i.e., one in each preconditioning step. The overall 
solution time is given as well. 

Results of the first experiment are listed in Table 1. It is important to clarify that 
beating the CSL or the CSG precondtioner is not the aim of this work. We rather focus on 
obtaining a preconditioner which can come in close comparison to them in performance, 
and is comparatively much easier to invert. The QD preconditioner takes around 3 times 
the number of iterations compared to the other choices. For this model problem (only) we 
also found that feeding in the preconditioner solution computed to a tolerance of 10~ 2 , 
to the Krylov method, as the starting guess, gives us a benefit of 50 iterations. 

The rest of the experiments are listed in Tables 2,3,4. They depict that the QD precon- 
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ditioner's performance comes within a factor of 3 of the other preconditioners even with 
strong spatial dependence in the wavenumber. The tuning effort is also much less, in 
fact, the relaxation parameters, the smoothing and the grid transfer methods can all stay 
constant. To date, the authors are not aware of any scientific method that minimizes the 
CSL shift or the CSG rotation angle for different problems, without extra overhead. Note 
that these tunable parameters have a pivotal role in establishing CSL/ CSG superiority in 
speed over the QD preconditioner. The experimental tables show the best cases for the 
CSL/ CSG preconditioners after they were hand-tuned for these parameters. This points 
to the possibility that the QD preconditioner might be used in an automatic solver setting. 
We tested the QD preconditioner against the Laplacian preconditioner (which can also 
be used in an automatic setting) and found the QD to be much superior in performance. 

For multidimensional Helmholtz operators (including the 2D operator), the critical 
eigenvalue A used in the QD preconditioner may be obtained from a one-dimensional 
counterpart, as is done in the current experiments. For Helmholtz problems with piece- 
wise constant wavenumber, the maximum discrete wavenumber value may be used as 
a rough approximation of Aq. However, this is also apt to bias the QD preconditioner 
spectrum more to the right than is really required. 



6 Conclusions and Outlook 



In this paper we showed that the Schrodinger equation for ionization problems can be 
decomposed into a coupled Helmholtz problem. This diagonal blocks of this coupled 
system consists of two-dimensional and three-dimensional Helmholtz problems. We 
propose Helmholtz model problems from these diagonal blocks. The blocks have ho- 
mogeneous dirichlet boundaries at one side and exterior complex scaling absorbing lay- 
ers (ECS-ABL) at the other side. Finite difference discretization (for non-uniform grids) 
results in a pitchfork-shaped spectrum which is largely distributed in the fourth quad- 
rant, but also has some parts crossing over in the third. Another property is that the 
spectrum is rather close to the real axis, and discrete problems are thus very challeng- 
ing to solve iteratively. We solved them iteratively using the preconditioned Bi-CGSTAB 
method and also presented the quadrant definite (QD) preconditioner, which we derive 
from a time integration scheme for the Schrodinger equation. We tried using GMRES 
and restarted GMRES but found that for the current problems these methods failed to 
reach their superlinear convergence phase. As a gross estimate we rate the efficiency of 
this preconditioner between the CSL/ CSG preconditioners and the Laplacian precondi- 
tioner, and it has the added advantage of having a multigrid favorable spectrum, i.e., its 
spectrum lies entirely in the fourth quadrant. This preconditioner can potentially be used 
in a automatic Helmholtz solver. The advantage of the QD preconditioner is that it can be 
built from standard multigrid components and it can be implemented matrix-free which 
significantly reduces the memory use. 
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Although we have used a low order discretization of the differential operators and a 
low order absorbing boundary conditions, we believe that calculations with higher order 
methods will lead to similar conclusions on the performance of the iterative method. 

Helmholtz problems, from an iterative perspective, can roughly be categorized into 
two classes which can be defined according to the available computational resources. 
One, where storing matrix operators is a possibility, and the other, where an iterative 
solution might have to be worked out using vectors alone. In the first situation, ILU(O) 
smoothing, and the Galerkin coarse grid operator used in a V(0,1) cycle render a very 
attractive multigrid method for preconditioner inversion. 

However, for the other class, the situation is comparatively much worse. First, for 
all preconditioners with a spectrum that leads to an efficient Krylov-subspace conver- 
gence, there is no appropriate smoother for multigrid. Second, we have to do with re- 
discretizing the Helmholtz operator on the coarse grid. This seems to work with non- 
standard F-cycles (with multiple coarse grid recursions), which are expensive. In future, 
we intend to investigate, how smoothing may be enhanced for matrix-free Helmholtz 
solution contexts, as well as how to bring multigrid down to work in V cycles for precon- 
ditioner inversion. 
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